{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "class HMM:\n",
    "    def __init__(self):\n",
    "        self.A=array()\n",
    "        self.B=array()\n",
    "        self.pi=array()\n",
    "        self.O=[] #观测序列\n",
    "        self.t=len(self.O)\n",
    "        self.m=len(self.A)\n",
    "        self.n=len(self.B)\n",
    "\n",
    "    def forward(self):\n",
    "        for i in range(self.m):\n",
    "            self.x[0][i]=self.pi[i]*self.B[i][self.O[0]]\n",
    "        for j in range(1,self.t):\n",
    "            for i in range(self.m):\n",
    "                temp=0\n",
    "                for k in range(self.m):\n",
    "                    temp+=self.x[j-1][k]*self.A[k][i]\n",
    "                self.x[j][i]=temp*self.B[i][self.O[j]]\n",
    "        result=0\n",
    "        for i in range(self.m):\n",
    "            result+=self.x[self.t-1][i]\n",
    "            print(\"前向概率矩阵及当前观测序列概率如下：\\n\")\n",
    "            print(self.x)\n",
    "            print(result)\n",
    "    # def get_stateprobability(self,t,p):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 使用第三方包实现HMM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 导入数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ID</th>\n",
       "      <th>NDVI</th>\n",
       "      <th>Forest</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.69</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>0.74</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>0.76</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>0.63</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>0.64</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1756</th>\n",
       "      <td>1757</td>\n",
       "      <td>0.17</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1757</th>\n",
       "      <td>1758</td>\n",
       "      <td>0.18</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1758</th>\n",
       "      <td>1759</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1759</th>\n",
       "      <td>1760</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1760</th>\n",
       "      <td>1761</td>\n",
       "      <td>0.13</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1761 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        ID  NDVI  Forest\n",
       "0        1  0.69       0\n",
       "1        2  0.74       1\n",
       "2        3  0.76       1\n",
       "3        4  0.63       0\n",
       "4        5  0.64       0\n",
       "...    ...   ...     ...\n",
       "1756  1757  0.17       0\n",
       "1757  1758  0.18       0\n",
       "1758  1759  0.16       0\n",
       "1759  1760  0.16       0\n",
       "1760  1761  0.13       0\n",
       "\n",
       "[1761 rows x 3 columns]"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from hmmlearn import hmm\n",
    "from collections import Counter\n",
    "from sklearn.model_selection import train_test_split,GridSearchCV\n",
    "from sklearn.metrics import accuracy_score,confusion_matrix\n",
    "file_path=\"D:\\Lenovo\\Desktop\\云南大学\\空间数据挖掘\\实验数据\\data9.csv\"\n",
    "data=pd.read_csv(file_path)\n",
    "observed_train, observed_test, hidden_train, hidden_test = train_test_split(\n",
    "    data['NDVI'], data['Forest'], stratify=data['Forest'],test_size=0.2, random_state=42)\n",
    "forest=data['Forest']\n",
    "ndvi=data['NDVI']\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 计算训练集的转移概率矩阵"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.79711452, 0.20288548],\n",
       "       [0.75503356, 0.24496644]])"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 定义NDVI离散化函数\n",
    "# 若不进行离散化，则会报错：ValueError: Symbol counts should be nonnegative integers\n",
    "# 原因：HMM模型要求每个状态的观测值数量必须是非负整数\n",
    "def discretize_ndvi(ndvi):\n",
    "    if ndvi < 0.1:\n",
    "        return 0\n",
    "    elif 0.1<=ndvi<0.3:\n",
    "        return 1\n",
    "    elif 0.3<=ndvi<0.5:\n",
    "        return 2\n",
    "    elif 0.5<=ndvi<0.7:\n",
    "        return 3\n",
    "    else:\n",
    "        return 4\n",
    "    \n",
    "observed_train_discrete = observed_train.apply(discretize_ndvi)\n",
    "observed_test_discrete = observed_test.apply(discretize_ndvi)\n",
    "\n",
    "# 初始化转移计数矩阵\n",
    "trans_counts = np.zeros((2, 2))\n",
    "\n",
    "# 计算状态转移频率\n",
    "for (i, j) in zip(hidden_train[:-1], hidden_train[1:]):\n",
    "    trans_counts[i, j] += 1\n",
    "\n",
    "# 转移概率矩阵\n",
    "trans_prob = trans_counts / trans_counts.sum(axis=1, keepdims=True)\n",
    "trans_prob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 计算观测状态概率矩阵"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.         0.14684685 0.2036036  0.57207207 0.07747748]\n",
      " [0.         0.         0.         0.27852349 0.72147651]]\n"
     ]
    }
   ],
   "source": [
    "# 初始化观测概率矩阵\n",
    "# 获取隐藏状态和观测状态\n",
    "hidden_states = hidden_train\n",
    "observations = observed_train_discrete\n",
    "n_hidden_states = 2 #隐藏状态有2，有森林或无森林\n",
    "n_observation_states = 5 #观测状态有5，NDVI重分类结果数\n",
    "emission_prob = np.zeros((n_hidden_states, n_observation_states))\n",
    "\n",
    "# 计算观测概率矩阵\n",
    "for hidden_state in range(n_hidden_states):\n",
    "    for observation_state in range(n_observation_states):\n",
    "        count = np.sum((hidden_states == hidden_state) & (observations == observation_state))\n",
    "        emission_prob[hidden_state, observation_state] = count\n",
    "\n",
    "# 将计数转换为概率\n",
    "emission_prob = emission_prob / emission_prob.sum(axis=1, keepdims=True)\n",
    "# emission_prob = (emission_prob + 1e-6) / (emission_prob.sum(axis=1, keepdims=True) + 1e-6 * n_observation_states)\n",
    "print(emission_prob)\n",
    "\n",
    "\n",
    "start_prob=np.array([1110/1408, 298/1408])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "MultinomialHMM has undergone major changes. The previous version was implementing a CategoricalHMM (a special case of MultinomialHMM). This new implementation follows the standard definition for a Multinomial distribution (e.g. as in https://en.wikipedia.org/wiki/Multinomial_distribution). See these issues for details:\n",
      "https://github.com/hmmlearn/hmmlearn/issues/335\n",
      "https://github.com/hmmlearn/hmmlearn/issues/340\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy_score:\n",
      "0.49008498583569404\n",
      "confusion_matrix:\n",
      "[[137 141]\n",
      " [ 39  36]]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# 定义HMM模型\n",
    "model=hmm.MultinomialHMM(n_iter=1000, n_components=2,init_params=\"est\",\n",
    "                         random_state=0)\n",
    "# 初始化参数\n",
    "# model.startprob_ = np.array([1110/1408, 298/1408]) #初始状态分布\n",
    "# model.transmat_ = trans_prob #隐藏状态转移矩阵\n",
    "# model.emissionprob_ = emission_prob #观测状态概率矩阵(发射矩阵)\n",
    "\n",
    "observed_train_discrete=np.array(observed_train_discrete)\n",
    "observed_test_discrete=np.array(observed_test_discrete)\n",
    "# fit函数只接受二维数组，reshape(-1, 1)用于将原有的一维数组转为n行1列的二维数组\n",
    "model.fit(observed_train_discrete.reshape(-1,1),len(observed_train))\n",
    "# 使用维特比算法解码最可能的隐藏状态序列\n",
    "hidden_states_pred=model.predict(observed_test_discrete.reshape(-1,1))\n",
    "log_probability=model.score(observed_test_discrete.reshape(-1,1))\n",
    "print(\"accuracy_score:\")\n",
    "print(accuracy_score(hidden_test, hidden_states_pred))\n",
    "print(\"confusion_matrix:\")\n",
    "print(confusion_matrix(hidden_test, hidden_states_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "n_iter的取值会显著影响模型精度，尝试使用网格搜索法找到最优n_iter，貌似网格搜索法不适用于HMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [],
   "source": [
    "# param_grid = {'n_iter': np.arange(10,100,1)}\n",
    "# # 使用 GridSearchCV 进行网格搜索\n",
    "# grid_model=hmm.MultinomialHMM(n_components=2, random_state=0, init_params='ste')\n",
    "# grid_search = GridSearchCV(grid_model, param_grid, cv=5,n_jobs=-1)\n",
    "# grid_search.fit(observed_train_discrete, hidden_train)\n",
    "# best_n_iter = grid_search.best_params_['n_iter']\n",
    "# best_score = grid_search.best_score_\n",
    "# print(best_n_iter)\n",
    "# print(best_score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 使用连续HMM测试"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7620396600566572\n",
      "[[230  48]\n",
      " [ 36  39]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\ProgramData\\anaconda3\\envs\\myenv2\\lib\\site-packages\\sklearn\\cluster\\_kmeans.py:1436: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=6.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "model2=hmm.GaussianHMM(n_iter=100, n_components=2,init_params=\"mc\",verbose=False,\n",
    "                       implementation=\"scaling\",params=\"mc\",random_state=0)\n",
    "model2.startprob_ = start_prob #初始状态分布\n",
    "model2.transmat_ = trans_prob #隐藏状态转移矩阵\n",
    "observed_train=np.array(observed_train)\n",
    "model2.fit(observed_train.reshape(-1,1))\n",
    "observed_test=np.array(observed_test)\n",
    "hidden_states_pred2=model2.predict(observed_test.reshape(-1,1))\n",
    "print(accuracy_score(hidden_test, hidden_states_pred2))\n",
    "hidden_states_pred2\n",
    "print(confusion_matrix(hidden_test,hidden_states_pred2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "myenv2",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
